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ABSTRACT 

Realistic models of magnetic reconnection in the solar chromosphere must take into 
account that the plasma is partially ionized and that plasma conditions within any 
two magnetic flux bundles undergoing reconnection may not be the same. Asymmetric 
reconnection in the chromosphere may occur when newly emerged flux interacts with 
pre-existing, overlying flux. We present 2.5D simulations of asymmetric reconnection 
in weakly ionized, reacting plasmas where the magnetic field strengths, ion and neutral 
densities, and temperatures are different in each upstream region. The plasma and 
neutral components are evolved separately to allow non-equilibrium ionization. As 
in previous simulations of chromospheric reconnection, the current sheet thins to the 
scale of the neutral-ion mean free path and the ion and neutral outflows are strongly 
coupled. However, the ion and neutral inflows are asymmetrically decoupled. In cases 
with magnetic asymmetry, a net flow of neutrals through the current sheet from the weak 
field (high density) upstream region into the strong field upstream region results from 
a neutral pressure gradient. Consequently, neutrals dragged along with the outflow 
are more likely to originate from the weak field region. The Hall effect leads to the 
development of a characteristic quadrupole magnetic field modified by asymmetry, but 
the X-point geometry expected during Hall reconnection does not occur. All simulations 
show the development of plasmoids after an initial laminar phase. 

Subject headings: magnetic reconnection — methods: numerical - plasmas — Sun: 
chromosphere 
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1. INTRODUCTION 


Magnetic reconnection in the lower solar atmosphere manifests itself through dynamical events 
such as jets (e.g., Shibata et al. 2007; Tian et al. 2014), explosive events such as Ellerman bombs 
(Ellerman 1917; Chen et al. 2001; Nelson et al. 2013; Yang et al. 2013; Peter et al. 2014), and possi¬ 
bly Type II spicules (De Pontieu et al. 2007). Reconnection events contribute to the heating of the 
chromosphere (e.g., Jess et al. 2014), and the contribution of jets and Type II spicules to the mass 
and energy budgets of the corona and solar wind is under active investigation (De Pontieu et al. 
2011; Cranmer &; van Ballegooijen 2010; Madjarska et al. 2011). Arge &; Mullan (1998) proposed 
that chromospheric reconnection could cause the elemental fractionation that is responsible for the 
first ionization potential (FIP) effect (see also Sturrock 1999; Feldman & Widing 2002). Partial 
ionization effects are important not just in the chromosphere but also in molecular clouds, pro¬ 
toplanetary disks, the neutral phases of the interstellar medium, some exoplanetary atmospheres 
(Koskinen et al. 2010), Earth’s ionosphere (Leake et al. 2014), the edge of tokamaks (Adams 2003), 
and dedicated reconnection experiments (e.g., Lawrence et al. 2013). 

Theories and simulations of reconnection have generally assumed that the plasma is fully 
ionized. This approximation is valid for the fully ionized solar corona but invalid for weakly ionized 
plasmas such as the solar chromosphere which has ionization fractions ranging from <0.01 to ^0.5. 
Partial ionization effects modify the dynamics of reconnection in several ways (Zweibel 1989; Zweibel 
et al. 2011). Before the onset of reconnection, current sheets thin significantly due to ambipolar 
diffusion (Brandenburg & Zweibel 1994, 1995). When there is strong coupling between ions and 
neutrals (e.g., on length scales longer than the neutral-ion mean free path, A n) ;), the effective ion 
mass is increased by the ratio of the total mass density p to the ion mass density pi. This decreases 
the bulk Alfven speed and consequently the predicted reconnection rate (Zweibel 1989). The plasma 
resistivity has contributions from both electron-ion and electron-neutral collisions (e.g., Piddington 
1954; Cowling 1956; Ni et al. 2007). 

The Hall effect is expected to be important on scales comparable to or less than the ion inertial 
length di . However, the effective ion inertial length is predicted to be enhanced in plasmas with 
strong ion-neutral coupling: d! i = di^/p/pi (e.g., Pandey & Wardle 2008). Consequently, it has 
been proposed that Hall reconnection may occur at longer length scales or larger densities than 
would be predicted from the ion density alone (Malyshkin & Zweibel 2011; Vekstein & Kusano 
2013). Malyshkin Sz Zweibel (2011) predict that enhancement of the Hall effect will lead to fast 
reconnection in molecular clouds and protoplanetary disks, but that this transition is less likely 
to be important in the solar chromosphere. In general, two regimes may be considered. When 
A ni << di, the Hall effect will be enhanced because of strong coupling and fast reconnection is 
predicted to occur. When d* < A ru < ri', some enhancement of the Hall effect due to ion-neutral 
coupling is also expected. However, as shown previously by Leake et al. (2012, hereafter, L12) and 
Leake et al. (2013, hereafter, L13), the ion density (and therefore also the local values for di, d! % , 
and X n i) may vary by more than an order of magnitude between the reconnection current sheet 
and the ambient plasma, which makes it difficult to a priori evaluate the importance of the Hall 
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effect in determining the rate of reconnection and the structure of the reconnection region in this 
parameter regime. 

The behavior of weakly ionized plasmas can be modeled using either single-fluid or multi-fluid 
formulations. A single-fluid approach incorporates the ambipolar diffusion term into the generalized 
Ohm’s law to account for relative drift between ions and neutrals (e.g., Brandenburg & Zweibel 
1994, 1995). This approach is less computationally expensive but implicitly assumes that the ions 
and neutrals are in ionization equilibrium. Multi-fluid models evolve each of the ion, neutral, and 
sometimes electron components separately to allow relative drifts between each of these populations 
and the fluid to be out of ionization equilibrium (e.g., Meier & Shumlak 2012). 

In recent years, considerable progress has been made in modeling magnetic reconnection in 
the lower solar atmosphere. Sakai et al. (2006), Smith & Sakai (2008), Sakai & Smith (2008), and 
Sakai & Smith (2009) performed two-fluid (ion-neutral) simulations of coalescing current loops in 
the solar chromosphere. Smith &; Sakai (2008) found that the rate of reconnection in simulations 
of the upper chromosphere was about ~20 times greater than in the lower chromosphere which 
has a considerably lower ionization fraction. Sakai & Smith (2008) and Sakai & Smith (2009) 
investigated the role of reconnection in penumbra filaments. These models assume fixed ionization 
and recombination rates rather than assuming a dependence on temperature and density. However, 
it is essential to incorporate rates that depend on physical conditions in order to accurately capture 
the important role of recombination during chromospheric reconnection. 

L12 and L13 used the plasma-neutral module of the HiFi framework (Lukin 2008; Meier 2011) 
to model chromospheric reconnection. The simulations showed a strong enhancement of ion density 
inside the current sheet as a result of ions being preferentially dragged along by the reconnecting 
magnetic field. Recombination and ion outflow were of comparable importance in the ion continuity 
equation for removing ions from the current sheet. The ions and neutral inflows were decoupled, 
but the ion and neutral outflow jets were strongly coupled. Late in time, the simulations by 
L12 and L13 showed the development of the secondary tearing instability known as the plasmoid 
instability (Loureiro et al. 2007). The increase in the ion and electron densities in the plasmoids led 
to an increase in the recombination rate that allowed further contraction of the magnetic islands 
on timescales comparable to the advection time of islands out of the current sheet. Ni et al. 
(2015) present one-fluid simulations of the plasmoid instability in the solar chromosphere using the 
NIRVANA code with and without guide fields. They find that the onset of the plasmoid instability 
leads to fast reconnection, and that slow shocks develop near the X-points. Cases without guide 
fields show rapid thinning of the current sheet due to ambipolar diffusion and radiative cooling, 
but these effects are not significant during guide field simulations. 

Prior models of magnetic reconnection in weakly ionized plasmas have generally assumed 
symmetric inflow. In general, however, there will be some asymmetry in the upstream magnetic 
field strengths, temperatures, and densities. The standard model for anenome jets predicts that 
chromospheric reconnection occurs when newly emerged flux interacts with pre-existing, overlying 
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flux (e.g., Shibata et al. 2007). It is reasonable to expect that these different plasma domains will 
have different plasma parameters, so this configuration naturally leads to asymmetric reconnection. 
Because the chromosphere is a dynamic magnetized environment (Leenaarts et al. 2007), asymmetry 
in the reconnection process is likely to be the norm. 

The physics of asymmetric inflow reconnection has been investigated in detail for fully ionized 
plasmas. One of the principal applications of this work has been Earth’s dayside magnetopause. 
Asymmetric inflow reconnection has also been investigated in the context of Earth’s magnetotail 
and elsewhere in the magnetosphere (Oieroset et al. 2004; Muzamil et al. 2014), the solar atmosphere 
(Nakamura et al. 2012; Murphy et al. 2012; Su & van Ballegooijen 2013; Su et al. 2013), laboratory 
experiments (Yamada 2007; Yoo et al. 2014; Murphy & Sovinec 2008), and plasma turbulence 
(Servidio et al. 2009, 2010). Cassak & Shay (2007) performed a scaling analysis for asymmetric 
inflow reconnection. They found that the outflow is governed by a hybrid upstream Alfven speed 
that is a function of the magnetic field strength and density in both upstream regions (see also 
Birn et al. 2010) and that the flow stagnation point in the simulation frame and magnetic field 
null were not colocated. The structure and dynamics of asymmetric reconnection in fully ionized 
collisionless and two-fluid plasmas have been studied previously in several works (e.g., Cassak & 
Shay 2008, 2009; Malakit et al. 2010, 2013; Pritchett 2008; Pritchett & Mozer 2009; Mozer et al. 
2008; Mozer &: Pritchett 2009; Swisdak et al. 2003; Aunai et al. 2013a,b; Murphy & Sovinec 2008). 
Simulations of the plasmoid instability during reconnection with asymmetric upstream magnetic 
fields by Murphy et al. (2013) showed that the resultant magnetic islands developed primarily into 
the weak field upstream region. Because the reconnection jets impacted the islands obliquely rather 
than directly, the islands developed net vorticity. In addition to asymmetric inflow reconnection, 
several groups have investigated asymmetric outflow reconnection (e.g., Oka et al. 2008; Murphy 
et al. 2010; Murphy 2010) and reconnection with three-dimensional asymmetry (e.g., Al-Hachami 
& Pontin 2010; Wyper & Jain 2013). 

Observational investigations of partially ionized reconnection in the chromosphere are chal¬ 
lenging because the dissipation length scales are significantly shorter than can be resolved with 
current instrumentation. Diagnosing the chromospheric magnetic field is an important but very 
difficult problem that requires inversion of spectropolarimetric data (e.g., Kleint 2012), and this will 
be even more challenging for short-lived dynamical events. Comparisons between simulations and 
observations should therefore concentrate on the large-scale consequences of chromospheric recon¬ 
nection which depend on small-scale processes. A complementary approach to solar observations 
is to study partially ionized reconnection in the laboratory. Lawrence et al. (2013) have performed 
such studies at the Magnetic Reconnection Experiment (MRX; Yamada et al. 1997). 

Our motivation is to investigate the role of asymmetry on the small-scale physics of reconnec¬ 
tion in partially ionized chromospheric plasmas. While there are similarities to symmetric cases, 
there are also physical effects such as neutral flows through the current sheet that do not occur in 
cases with symmetric inflow. We use the plasma-neutral module of the HiFi modeling framework 
(Lukin 2008) to perform simulations of magnetic reconnection with asymmetric upstream magnetic 
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field strengths, temperatures, and/or densities in the weakly ionized solar chromosphere. In Section 
2, we describe the numerical method and problem setup (see also L12 and L13). In Section 3, we 
describe the global dynamics of reconnection, the structure of the current sheet, the role of the Hall 
effect, the dynamics of the plasmoid instability, the motion of the X-point, and ion/neutral flows 
at the X-point. We discuss our results in Section 4. 


2. NUMERICAL MODEL AND PROBLEM SETUP 

The HiFi framework 2 (Lukin 2008; Glasser Sz Tang 2004) uses a spectral element spatial 
representation and implicit time advance to solve systems of partial differential equations. The 
modular approach makes new physical models straightforward to implement (e.g., Gray et al. 2010; 
Lukin & Linton 2011; Ohia et al. 2012; Le et al. 2013; Stanier et al. 2013; Browning et al. 2014; Lee 
et al. 2014). In this paper, we use the module for partially ionized, reacting plasmas (Meier 2011, 
L12, L13). These plasmas consist of neutral and ionized hydrogen and electrons. A multi-fluid 
approach allows the plasma and neutral components to be modeled separately. We summarize the 
equations used for our 2.5D simulations in this section, but direct the reader to L12 and L13 for 
further detail. In general, we use the notation and conventions from L12 and L13. The subscripts 
‘n’, ‘i’, and ‘e’ refer to neutrals, ions, and electrons, respectively. 


2.1. Normalizations 

Following L12 and L13, we normalize the fluid equations to values characteristic of the lower 
chromosphere. We choose a characteristic length scale of L* = 1 x 10 4 m, a characteristic number 
density of n* = 3 x 10 16 m -3 , and a characteristic magnetic field strength of B* = 1 x 10~ 3 T. From 
these quantities, we derive additional normalizing values to be U* = B*/ m p n* = 1.26 x 10 5 
ms -1 for velocity, f* = L*/U* = 0.0794 s for time, T* = B 2 //cb/Uo?U = l-92x 10 6 K for temperature, 
R* = B 2 //i 0 = 0.796 Pa for pressure, J* = B*//ioR* = 7.96 x 10” 2 A-m~ 2 for current density, and 
V* = HoL*V± = 1.58 x 10 3 0 • m for resistivity. Unless otherwise indicated (e.g., Sections 2.2 and 
2.4), the equations presented for the simulation will be in dimensionless units according to these 
normalizations. 


2.2. Ionization and Recombination 

The HiFi module for partially ionized plasmas includes both ionization and recombination of 
hydrogen to allow departures from ionization equilibrium (L12; L13). The ionization and recombi- 


2 See http://f acuity.washingt on.edu/vlukin/HiFi_Framework.html 
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nation rates are given by 


rr - —n n v i ' 

r[ ec = -mv rt 


(1) 

( 2 ) 


with r* on = — F l ° n and Tt ec = —r^ ec . The ionization frequency of hydrogen is approximated to be 


n e A 


v = 


X + (frion/T* 


'J n * 


K 


exp 


0 / 


m 3 s _1 , 


( 3 ) 


with A = 2.91 x 10 -14 , K = 0.39, X = 0.232, and 0j on = 13.6 eV (Voronov 1997). Here, T* 
is the electron temperature in eV. Smirnov (2003) approximates the recombination frequency of 
hydrogen to be 

v rec = 2.6 x 10- 19 ^= m 3 s' 1 . (4) 

Vt! 

At a characteristic temperature of 9000 K, this corresponds to an equilibrium ionization fraction 
of 4.1 x 10“ 4 . These expressions do not include the consequences of non-local thermodynamic 
equilibrium radiative transfer. For example, if Lyman a is optically thick, then there will be more 
neutral hydrogen in the n = 2 state than predicted from these expressions which would allow for 
enhanced ionization. At lower temperatures, ionization from low-FIP elements may contribute 
more to the electron density than hydrogen. 


2.3. Multi-fluid Equations 

In this section, we summarize the equations solved by the plasma-neutral module of HiFi. The 
equations are described more thoroughly by L12 and L13 (see also Meier 2011; Meier & Shumlak 
2012). Our simulations include the modifications and extensions to the model described by L13. 


The ion and neutral continuity equations are 



^ + V-(niVi) = 

-prec _j_ -pion 

( 5 ) 

^ + V'Kv n ) = 

prec _j_ -pion 
n ' n ’ 

(6) 


These equations include ionization and radiative recombination as described in Section 2.2, L12, 
and L13. We assume quasineutrality such that m = n e . 

The ion and neutral momentum equations are given by 
<9 

— ( uiiriiVi ) + V • (///,//, V;V; + P* + P e ) = 

J x B + Rj" + T^rrnVn - T r n ec m. i y i + T cx mi (V n - V,) + - R“, (7) 

( rai ""V») + v ' (i7!jn„V„V„ + P„) = 

-Rf + r r n ec miVi - r i on m t V n + T cx mi (Vi - V n ) - R^ + R“. 


( 8 ) 
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The Lorentz force acts directly on the plasma but not the neutrals, while the neutral pressure 
gradient acts directly on the neutrals but not the plasma. Coupling between the plasma and 
neutrals is achieved though many of the terms on the right hand side of Eqs. 7 and 8. The term 
Rm re p resen ts momentum transfer from neutrals to ions due to identity preserving collisions and 
is given by 

RT = m in niv in (V n - Vi ), (9) 

where = mim n /(mi+m n ) and the collision frequency I'jn is given by Eq. 7 of L13. The pressure 
tensor for species a is 

P a = P a I + 7T q , (10) 

where P a is the scalar pressure and I is the identity tensor. The viscous stress tensor is then 

vr Q = [W a + (VV a ) T ] , (11) 

with as the isotropic dynamic viscosity coefficient. The terms that depend on the ionization and 
recombination rates represent momentum transfer by particles that change identity from neutrals to 
ions and vice-versa. We follow L13 and include charge exchange. Here, T cx is the charge exchange 
reaction rate and R^ is the momentum transfer due to charge exchange rections from species /3 
to species a. Charge exchange leads to increased momentum and energy transfer between species 
by about a factor of two, and therefore more effective coupling. We use the formulation for charge 
exchange given by Eqs. 4-6 of L13 (see also Meier 2011; Meier & Shumlak 2012; Leake et al. 2012; 
Barnett et al. 1990). 

We adjust the collisional cross sections to be E m = = 5 x 10 -19 m 2 (used in Eq. 7 of L13 to 

calculate Uj /a ) and E nn = 5 x ICE 19 m 2 (used in Eq. 9 of L13 to calculate v nn )- These values differ 
from L12 and L13, but the value for is consistent with Khomenko & Collados (2012) and Ni 
et al. (2015). These cross sections are a function of energy, but we use constant values appropriate 
for the chromosphere as a simplifying assumption (see also Draine et al. 1983). The neutral and ion 
viscosity coefficients £ n and are set by Eq. 8 of L13 using the revised value of E nn . While is a 
function of local physical conditions, and £ e are constant throughout the domain and correspond 
to the parallel component of the Braginskii viscosity for each species computed using the mean of 
the asymptotic upstream densities and temperatures. We use = 2.8 x 10 -6 and £ e = 4.6 x ICE 14 
with a normalization of £* = m p n*L* 14 - 

The neutral-ion mean free path 

A ni = ^ (12) 

V • 

ni 

governs the lengths scales above which the neutrals and ions are coupled to each other. Here, 
the neutral thermal velocity is Vr, n = \/2 kBT n /m n where ks is the Boltzmann constant, the 
neutral temperature is T n , and the neutral mass is m n . The frequency v' ni = v n i + includes 
contributions from the neutral-ion collision frequency lyu (see Eq. 16 of L12) and the neutral-ion 
charge exchange frequency i/W (see also Eqs. 4-6 of L13). The ions and neutrals will be coupled 
on length scales much longer than A m , and decoupled on scales much shorter than A m . 
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The energy equations for the plasma and neutral components are given by Eqs. 19 and 20 
of L12. Both equations include frictional heating due to identity preserving collisions, thermal 
transfer due to changes in identity, and the effects of charge exchange. The energy equation for the 
plasma component combines the electron and ion energy equations and includes Ohmic heating, 
optically thin radiative losses, and anisotropic thermal conduction (Eq. 21 of L12). The neutral 
energy equation includes isotropic thermal conduction that depends on plasma parameters and is 
set by Eq. 10 of L13. Neutral thermal conduction dominates thermal diffusion, in part due to 
rapid thermal transfer between neutrals and ions. For characteristic values of T = 9000 K and 
n n = 7.6 x 10 18 m -3 , the neutral thermal conductivity is n n = 3.05 x 10 22 m _1 s _1 . We assume 
that the ion and electron temperatures are equal, but the neutral temperature is evolved separately. 


The generalized Ohm’s law for this paper is given by 

„ , T , JxB VP e m e u en 

E + V i xB = r/J+ e e en 


(Vi - V n 


en e en e e 

The resistivity includes both electron-ion and electron-neutral collisions and is given by 


(13) 


n = 


m. e n e (u ei + Ven) 
(en e ) 2 


(14) 


where the electron-ion collision frequency v e i and the electron-neutral collision frequency v en are 
functions of number density and temperature and are given by Eq. 13 of L13 with E en = T, ne = 
1 x 10~ 19 m 2 as used previously in L12, L13, Khomenko & Collados (2012), and Ni et al. (2015). 
The resistivity therefore depends on plasma parameters and is not a constant as in L12. We include 
the Hall term and evolve scalar plasma and neutral pressures. 


2.4. Initial Conditions 


Here we describe the procedure to establish an approximate initial equilibrium with asym¬ 
metric upstream magnetic field strengths, densities, temperatures, and ionization fractions. Before 
proceeding, we define several parameters to quantify asymmetries in different fields. The magnetic, 
temperature, number density (of ions and neutrals), and ionization fraction asymmetries are defined 
as 


B 

T 

M 


B2 
B\ ’ 

Ti 

T\ ’ 


+ T^n, 2 
1^1,1 + ^n,l 

h 

h' 


(15) 

(16) 
(17) 




( 18 ) 
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where the subscripts ‘ 1 ’ and ‘ 2 ’ correspond to the asymptotic upstream magnitudes of each field 
for y > 0 and y < 0 , respectively, and the ionization fraction is defined as 

rii 


f = 


(19) 


rii + n n 

All of these quantities are functions of time, so the subscript ‘0’ in expressions below indicates 
correspondence to the initial conditions. 


The equilibrium magnetic field is specified as a modified Harris sheet, 

B x o (y) = -E>i,o 


tanh ( — 6 ) + 6 


1 + 6 


( 20 ) 


where A^, is the initial thickness of the current sheet (see also Birn et al. 2008, 2010; Murphy et al. 
2012, 2013). The initial magnetic asymmetry is given by Bq = (1 — b)/(1 + 6 ) with the convention 
that 0 < b < 1 so that < Bgo- We describe the in-plane magnetic field using the magnetic 
flux, A z , such that 

B = V x (A z z) + B z z. (21) 

The flux corresponding to Eq. 20 is 

-£>i,o 


A z o(y) = 


1 + 6 


Av, In cosh ( — b ) + by 

Alp 


( 22 ) 


The out-of-plane current density is then 


J z o(y) = ~ 


B i,o 

/+A xp 


sech 2 ( +— 6 
1 + 6 


(23) 


The initial temperature is set using the relations 


T 0 (y) = T 1;0 [1 + (To - 1) (1-C 2 )], 


C = 2 


1 + tanh ( -^— 6 


(24) 

(25) 


The initial ionization fraction is found numerically by equating the ionization rate T“ n with the 
recombination rate Y[ ec using Eqs. 1-4. The initial conditions are therefore in ionization equilib¬ 
rium. If the initial temperature is non-uniform, then the initial ionization fraction will also be 
non-uniform: if 7o 7 ^ 1, then J -0 7 ^ 1. This is in contrast to L12 and L13 which both assume that 
the temperature and ionization fraction are both initially uniform. 

We define P tot to be the sum of the plasma and neutral pressures and the plasma pressure P p 
to be the sum of the electron an ion pressures, 


-Ptot — Bn + Bp, 
Pp = Pi + Be- 


(26) 

(27) 
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For equal temperatures, the neutral and plasma pressures are related to the ionization fraction by 


Pn 

P P 


Ptot 

Ptot 



(28) 

(29) 


where we recall that the total number of particles depends on the ionization fraction and that 
electrons are included in the plasma pressure. We calculate P t ot,o to balance the Lorentz force 
associated with the magnetic field profile given by Eq. 20, 


Ptot,o(v) 


(1 + /5i,o) 


Bio 

2 /i 0 


B x o(y) 2 

2fi 0 


(30) 


where /3po is the ratio of Ttot,o to the magnetic pressure for y > 0. The number densities are then 
given by 


“ k B T' 

P p 

Ul 2k B T' 

We assume that the ions and electrons have equal temperatures and number densities. 


(31) 

(32) 


The above magnetic and pressure profiles satisfy the relation 


0 — —V ( Pp + P n ) + J x B. (33) 

However, there must also be a relative velocity between the plasma and neutrals to allow coupling 
via the momentum equations of the different species through collisions and charge exchange. In 
the absence of charge exchange, the steady state momentum equations for the plasma and neutral 
components are 


0 = -VP p + J x B + Rf\ (34) 

0 = ~VP n - Rf\ (35) 


respectively, where the frictional force R^ n is defined in Eq. 9. In this approximation, the initial 
ion and neutral velocities along the inflow direction can be chosen so that the initial configuration 
contains no net force on either of the plasma and neutral components, 


WyO PnyO — 


dPn/dy 


where the neutral pressure gradient is calculated analytically using the expression 

dP n (1 — fo\ ( Bi,o \ n / \_i 2 ( V 


dy 


l + /o 


, 0 x*li + b )) BMsech: \x, 


-b 


(36) 


(37) 
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Case E: Initial Conditions 



Fig. 1.— The initial conditions for case E with Bq = 1, 7o = 0.95, Mq = 1.47, and J-q = 0.35. 

Because our simulations include charge exchange, momentum transfer between the ions and neutrals 
is of order twice as effective as the case with frictional coupling alone. In lieu of an exact solution 
to the steady state momentum equation, we reduce the initial ion-neutral drift velocity given in Eq. 
36 by a factor of two so that the initial conditions for both ions and neutrals are in approximate but 
not exact force balance. We set Vi y o and V ny o so that they have the same magnitude but opposite 
sign. Example initial conditions are shown in Fig. 1. 

The initial state does not represent an exact equilibrium, which leads to an outward moving 
pulse along the inflow direction that propagates from the initial current sheet. This pulse is mostly 
damped by placing a layer of significantly enhanced viscosity near the conducting walls at y = ±1. 
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2.5. Boundary Conditions 

We simulate a half-domain that extends from 0 < x < L x and from —L y < y < L y , where x is 
the outflow direction and y is the inflow direction. We assume that the fields A z , rii, n n , pi, p n , Vi y , 
V ny , Vi z , V nz , and J z are symmetric about x = 0 [e.g., A z (x,y) = A z (—x,y)], and that the fields 
V x and B~ are antisymmetric about x = 0 [e.g., V x (x,y) = —V x (—x,y)\. We assume that all fields 
are periodic along the x direction with a period of 2 L x . We apply perfectly conducting, zero-flux 
boundary conditions along the boundaries at y = ±L y . 

The choice of initial and boundary conditions enforces that the reconnection outflow is sym¬ 
metric about x = 0. The development of the plasmoid instability under these conditions often leads 
to the formation of a magnetic island centered about x = 0 that cannot be advected out of the 
system. In general, there will be some degree of asymmetry along the outflow direction that can 
modify the internal structure of the reconnection layer and the dynamics of reconnection (Oka et al. 
2008; Murphy et al. 2010; Murphy 2010). This can be achieved by simulating the whole domain of 
interest and allowing the driving term, the outflow boundaries, or the initial perturbations to be 
asymmetric (e.g., Murphy et al. 2013). 


2.6. Initializing Reconnection 


To initialize reconnection, we apply a source function to the evolution equation for magnetic 
flux of the form 

d A 

—^ = eA</,A (x, y) Q (x, y) 0 (t) (38) 

where 


A (x,y) 

tt(x,y) 

Q(t) 



(39) 

(40) 

(41) 


The spatial length scales for this source function are given by h x = 4A^ and h y = A^. The 
function A(x, y) localizes the source function in a form akin to a tearing mode eigenfunction to 
create an X-point at the field reversal along x = 0. The function fl(x, y) ensures that the pulse 
goes to zero along the outer boundaries. The pulse is applied between 0 < t < to with a waveform 
given by 0(t). We use e = 0.05 and to = 5. The application of this electric field ceases before 
the reconnection layer has had a chance to develop. In contrast, L12 applied a small, localized 
perturbation to the magnetic flux while L13 applied a small amplitude localized rotational flow 
perturbation to initialize reconnection. 
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3. LOCAL AND GLOBAL DYNAMICS OF RECONNECTION 

We present five simulations to investigate the impact of asymmetry during magnetic recon¬ 
nection in partially ionized chromospheric plasmas. The simulation parameters are shown in Ta¬ 
ble 1. Case A is the symmetric test case. Case B has asymmetric upstream temperatures, and 
consequently asymmetric upstream densities and ionization fractions, but symmetric upstream 
magnetic field strengths. Case C has symmetric upstream temperatures but a factor of two 
difference in the upstream magnetic field strengths. Cases D and E have asymmetric tempera¬ 
tures and magnetic field strengths. We keep three parameters constant between runs: the sum 
of the initial magnetic pressure, neutral pressure, and plasma pressure which equals 1.495 (1.19 
Pa); the mean of the initial asymptotic upstream magnetic energy densities which is given by 
\ (B'{ 0 /2 + B\ J 2) = 0.3125 (0.249 J m 3 ); and the mean of the initial asymptotic upstream 
temperatures which is (T lj0 + T2 i o)/2 = 4.69 x 10~ 3 (9000 K). 

The domain size for all simulations is ( L x , L y ) = (2,1). The resolution in Case A is rn x = 256 
elements along the outflow direction and rri y = 128 elements along the inflow direction. The 
resolution in Cases B-E is rn x = m y = 256. We use sixth order basis functions for all simulations, 
resulting in effective total resolution of ( M x , M y ) = 6(m x ,m v ). Grid packing is used to concentrate 
mesh in the reconnection region. In Case A, the current sheet does not move from y = 0 so the 
mesh packing along the inflow direction is concentrated to a thin region near y = 0. In Cases B D, 
the current sheet drifts slowly in the —y direction so the highest resolution is needed over a much 
longer distance to resolve the dynamics (see also Murphy et al. 2012, 2013). High resolution along 
the outflow direction helps capture the front end of the reconnection jet as well as the dynamics 
of the plasmoid instability. The resolution along the outflow direction from 0 < x < 1.2 and 
1.95 < x < 2 is approximately twice as high as from 1.2 < x < 1.95. 


3.1. Structure of Reconnection Region 

All five simulations show broadly similar evolution. The electric field application from t = 0 
to t = 5 described in Section 2.6 allows the inflow/outflow pattern associated with two-dimensional 
reconnection to develop in the portion of the current sheet near the origin. The outflow jet lengthens 
as it plows into a region of enhanced ion density associated with current sheet thinning outside of 
the reconnection region. By t ~ 25, laminar reconnection is well-established. The current sheet has 
thinned from its initial thickness of 0.1 to a thickness of 5 ~ 1.0 x 10 -3 (~10 m). This is comparable 
to the value of the neutral-ion mean free path evaluated inside the current sheet: A n j > 6.0 x 10~ 4 
to 1.0 x 10 ” 3 (>6 to 10 m). The ionization fraction inside each current sheet is of order 0.01 
at this time. The structure of the reconnection region for all simulations at t = 25 is shown in 
Figures 2, 3, 4, and 5. At t = 25, the current sheets have thinned to ~A n ,, laminar reconnection is 
well-established, and plasmoid formation has not yet begun. The plasmoid instability onsets during 
all of these simulations. The simulations end when structures develop on scales comparable to the 



Table 1. Simulation initial parameters 


Case 

-E>i,o 

-£>2,0 

Ti, 0 

T‘2,0 

m,o 

n 2 ,0 

/l,0 

f'2,0 

Pifi 

@2,0 


^ni, 2,0 

A 

7.9 

7.9 

9000 

9000 

7.6 

7.6 

0.00041 

0.00041 

3.8 

3.8 

220 

220 

B 

7.9 

7.9 

8750 

9250 

7.8 

7.4 

0.00024 

0.00068 

3.8 

3.8 

366 

136 

C 

10 

5 

9000 

9000 

00 

bo 

6.4 

0.00041 

0.00041 

11.0 

2.0 

261 

190 

D 

10 

5 

8750 

9250 

6.6 

8.6 

0.00024 

0.00068 

2.0 

11.0 

432 

117 

E 

10 

5 

9250 

8750 

6.2 

9.1 

0.00068 

0.00024 

2.0 

11.0 

163 

314 


Note. The units are G for magnetic field, K for temperature, 10 18 m -3 for number density 
(including both neutrals and ions), and m for the neutral-ion mean free paths. The neutral-ion 
mean free path A n j is calculated using Eq. 12 and includes charge exchange reactions. When 
calculating the charge exchange cross section using Eqs. 5-6 of L13, we use that the relative 
velocity between ions and neutrals is much less than the neutral and ion thermal speeds for the 
initial conditions. 
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resolution scale as a result of the plasmoid instability. 

Figure 2 shows the out-of-plane current density and the ion outflow. In Case B, the current 
sheet is slightly arched so that the X-point is closer to the low temperature upstream region (y > 0) 
than the ends of the current sheet. In the cases with magnetic asymmetry, the current sheet is 
arched so that the X-point is closer to the strong field upstream region than the ends of the outflow 
jets (see also Fig. 5b). Case E, which has the strong field side coincident with the high temperature 
side, shows the fastest development of these three cases and the strongest current density. 

As in L12 and L13, the ion and neutral outflows are tightly coupled. The right half of each 
panel in Figure 2 only show pseudocolor maps of V lx because V nx is very similar in the outflow 
jet. The strong coupling of the ion and neutral outflow jets occurs because the mean free path 
for neutral-ion collisions is comparable to the current sheet thickness which is significantly shorter 
than the length of the current sheet. 

Figure 3 shows the inflow components of the ion velocity on the left side of each panel and 
the neutral inflow speed on the right side (see also Figs. 5d and 5e). As in L12 and L13, the 
inflow velocities are decoupled. In Cases B -E, the decoupling is asymmetric in part because \ n { is 
different in each upstream region. The neutrals and the ions at the X-point are moving in opposite 
directions. For simulations with magnetic asymmetry (Cases C, D, and E), there is neutral flow 
through the current sheet from the weak magnetic field side toward the strong magnetic field side. 
This corresponds to a neutral pressure gradient that is pushing neutrals from the weak field side 
into the strong field side (see Fig. 5f). 

The ion density is strongly peaked within the current sheet for all cases, as shown on the left 
side of each panel in Fig. 4 (see also L12; L13). The decoupling of ions and neutrals on scales below 
the neutral-ion mean free path allows the ions to be swept into the current sheet by the magnetic 
field. Neutrals are dragged along by collisions with ions, though less effectively on these short length 
scales. The current sheet is therefore out of ionization equilibrium: the ionization fraction is much 
higher than the equilibrium value. As a result, recombination becomes of comparable importance 
to the outflow in the ion continuity equation. 

Next we consider the consequences of changing the temperature asymmetry while maintaining 
the same magnetic asymmetry. As we go from D to C to E, the magnetic asymmetry remains 
constant (B = 0.5) but the temperature asymmetry goes from T = 1.06 to 1 to 0.95. Case D has 
higher temperature in the weak field upstream region; Case C has initially uniform temperature; 
and Case E has higher temperature in the strong field upstream region. 

Switching the temperature asymmetry while maintaining the same magnetic asymmetry has 
a significant impact on how quickly reconnection develops, the reconnection rate at a given time, 
and the onset time and mode structure of the plasmoid instability. Reconnection develops more 
quickly and occurs at a faster rate in Case E where the strong field upstream region has a higher 
temperature than the weak field upstream region. As we proceed from Case D to C to E, the 
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(a) Case A at t=25 
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Fig. 2.— The out-of-plane current density J z (left) and the ion outflow (right) at t = 25 for Cases 
A-E. The right half of each plot shows color contours of the outflow component of ion velocity V zx , 
ion velocity vectors, and contours of magnetic flux A z . This image is scaled significantly along the 
y direction which exaggerates the inflow y-component of velocity. 



(e) Case E at t=25 



(d) Case D at t=25 




(b) Case B at t=25 
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(e) Case E at t=25 



(d) Case D at t=25 



(c) Case C at t=25 



(b) Case B at t=25 



(a) Case A at t=25 
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Fig. 3.— The inflow components of the ion velocity Vi y (left) and the neutral velocity V ny (right) in 
the reconnection region at t = 25 for Cases A-E. The solid black contours represent the magnetic 
flux A z . The dashed green contour indicates the locations where the inflow component of velocity 
equals zero. 
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Fig. 4.— The logarithm of the ion number density log 10 rii (left) and the out-of-plane magnetic 
field B z (right) at t = 25 for Cases A-E. The solid black contours represent the magnetic flux A z . 
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(e) Case E at t=25 



(d) Case D at t=25 
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y-yn 


Fig. 5.— The structure of the current sheet along x = 0 at t = 25 very near the current sheet 
relative to the X-point. Shown are (a) the reconnecting component of the magnetic field B x , (b) the 
out-of-plane current density J z , (c) the logarithm of ion density log 10 rn, (d) the inflow component 
of ion velocity Vi y , (e) the inflow component of neutral velocity V nyi and (f) the neutral pressure. 
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y-Vn 


Fig. 6.— The logarithm of the neutral-ion mean free path A n i = Vr.n/^li along x = 0 at t = 25. 
The minimum values of A n i in this plot are 6.5 x 10 -4 , 6.3 x 10 -4 , 7.9 x 10~ 4 , 7.2 x 10 -4 , and 
8.3 x 10” 4 for cases A through E, respectively. 

relative velocity between the ions and neutrals increases in magnitude on the weak field side but 
does not change substantially on the strong field side (see, e.g., Fig. 5d while noting that the V ny 
profiles in Fig. 5e are nearly identical between these three runs). The neutral-ion mean free path 
A n i along x = 0 at t = 25 is shown in Fig. 6. In this figure, the minimum value of A n i increases by 
10% from Case D to C and another 5% from Case C to E. In contrast, the values of A n i increase by 
~30% to ~60% from D to C to E in both near-upstream regions. The reconnection rate at t = 25 
measured as dA~/df evaluated at the X-point along x = 0 is 3.4 x 10~ 4 , 3.8 x 10~ 4 , and 4.5 x 10 -4 in 
Cases D, C, and E, respectively, which suggests that differences in ion-neutral coupling inside and 
around the current sheet are largely responsible for the observed differences between simulations. 
In particular, the reconnection rate in these simulations is greater when there is weaker coupling 
between neutrals and ions. The V ny profile in Fig. 5d remains nearly unchanged between simulations 
with the same magnetic asymmetry, which suggests that the neutral flow across the current sheet 
depends strongly on magnetic asymmetry but only weakly on temperature asymmetry. 


3.2. Role of the Hall Effect 

The Hall effect is not expected to be significant in the solar chromosphere during magnetic 
reconnection unless structures develop on scales comparable to or below the ion inertial length di 
(Malyshkin Sz Zweibel 2011). In the case of strong coupling between ions and neutrals, it has been 
proposed that the ion inertial length may be enhanced because each ion will be dragging along a 
much larger mass (Pandey & Wardle 2008; Malyshkin & Zweibel 2011); consequently the effective 
ion inertial length is proposed to be d[ = di^p/pi. We therefore include the Hall effect in these 
simulations. 
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The characteristic thickness of the current sheet is ~10 -3 (10 m; see Fig. 5). The ion inertial 
length in the upstream regions is initially of order ci; ~ 4 x ICC 4 (4 m; using upstream parameters 
from Case A), while d[ ~ 0.02 (200 m). In the current sheet, the pileup of ions leads to a local ion 
inertial length of di ~ 8 x ICC 5 (0.8 m; based on an ion density of ~3no), which corresponds to 
d[ ~ 8 x ICC 4 (8 m). These lengths scales can also be compared to \ ni which is of order ~0.01 to 
0.03 in the regions just upstream of the current sheet and ~1CC 3 in the current sheet. 

The right side of Figure 4 shows the quadrupole structure of the out-of-plane magnetic field 
B z that forms as a result of the Hall effect during laminar reconnection. The quadrupole field 
strength remains about an order of magnitude smaller than the asymptotic values of B x , but 
\B Z \/B x + By locally reaches ~0.5 near the location where \B Z \ is largest on the weak field 
upstream region in Cases C-E. This occurs in part because the in-plane field is weaker in the near 
upstream regions than in the far upstream regions. Despite the locally high value of the ratio of 
\B~\ to the in-plane field, these cases show neither the significant enhancement of the reconnection 
rate nor the development of the low aspect ratio (X-point) geometry expected during Hall-mediated 
reconnection in fully ionized plasmas (e.g., Biskamp et al. 1997). 

The structure of the quadrupole field is modified by both temperature and magnetic asym¬ 
metries. When the temperature or magnetic asymmetries are changed, this also impacts the neu¬ 
tral density, ion density, and ionization fraction asymmetries and therefore A n j. In Case B, the 
quadrupole lobe on the low ion density side (y < 0) has a greater spatial extent and higher mag¬ 
nitude than the high ion density side. In Cases C-E with magnetic asymmetry, the quadrupole 
lobe on the strong magnetic field side is localized near the current sheet while the quadrupole lobe 
on the weak magnetic field side has a greater spatial extent. This asymmetry of the quadrupole 
field is qualitatively similar with previous fully kinetic simulations (see, e.g., Fig. 8c. of Pritchett 
2008). As we proceed from Case D to C to E at t = 25, the quadrupole field on the strong field 
side increases in strength while the quadrupole field on the weak field side decreases in magnitude 
but increases in spatial extent. 

Even though the ion inertial scale is shorter than the current sheet thickness throughout these 
simulations, this might not be the case during the long-term evolution of the plasmoid instability. 
While simulations of fully ionized plasmas have shown that the plasmoid instability allows recon¬ 
nection to occur at a rate that is roughly independent of Lundquist number, Daughton et al. (2009) 
and Shepherd & Cassak (2010) have proposed that the most important role of the plasmoid insta¬ 
bility might actually be to allow structure to develop on scales smaller than the ion inertial length 
which would then allow fast collisionless reconnection to take place. This proposed mechanism 
has not yet been tested in weakly ionized plasmas, which would require capturing highly nonlinear 
behavior at later times and will be explored in future work. 
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Fig. 7.— The ion density rii (left) and current density J z (right) for Cases A-E near the end of 
each simulation after plasmoids form. 
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3.3. Plasmoid Formation 

Each simulation shows the formation of plasmoids after the current sheets have thinned suf¬ 
ficiently. Figure 7 shows the ion density, out-of-plane current density, and magnetic flux during 
the early nonlinear evolution of the plasmoid instability for each case. In contrast to Figures 2-6, 
the times for each panel were chosen to be near the end of each simulation but before the plasma 
substructures associated with the cascading nonlinear evolution of the plasmoids reached the grid 
scale. The appearance of additional in-plane null points occurs at t = 27.9 for Case A, t = 29.8 
for Case B, t = 29.6 for Case C, t = 31.1 for Case D, and t = 25.9 for Case E. In all cases, the 
secondary islands develop in the central portion of the current sheet: within x < 0.16 compared to 
a current sheet length of ~0.8. The out-of-plane current density has a local extremum near each 
X-point. 

We first consider the effect of introducing a temperature asymmetry by comparing Cases A and 
B. The plasmoid instability takes longer to develop in Case B which has temperature asymmetry, 
but the mode structure does not change substantially. Both cases show a central O-point with two 
X-points on each side with an additional X-point/O-point pair further out. The symmetry about 
x = 0 allows a pitchfork bifurcation to change the central X-point into an O-point. The resulting 
island then grows due to reconnection and is unable to be advected out of the reconnection region 
by the outflow. Prior simulations of the plasmoid instability in resistive MHD have avoided this 
situation by simulating the entire domain and including a slight asymmetry along the outflow 
direction so that any central islands that form do not remain in the current sheet indefinitely. 

The structure of the current sheet and the resultant plasmoid instability are more directly 
modified by magnetic asymmetry, including when it is coupled with temperature asymmetry. While 
Case D also undergoes a pitchfork bifurcation to change the central X-point into an O-point, Cases 
C and E maintain a central X-point and develop multiple alternating X-points and magnetic islands. 
As in prior simulations using the resistive MHD approximation (Murphy et al. 2013), the resulting 
islands preferentially grow into the weak field upstream region. 

We now compare our simulations with the complementary simulations of the plasmoid insta¬ 
bility during partially ionized chromospheric reconnection by Ni et al. (2015). Ni et al. capture 
the long term nonlinear evolution of the plasmoid instability using a single fluid approach that 
incorporates ambipolar diffusion and assumes ionization equilibrium. The parameter regimes for 
the two sets of simulations differ. For example, Ni et al. assume an initial upstream /3 of 0.1 which 
allows heating up to ~8 x 10 4 K, while our simulations have high f3 and thus do not result in 
significant heating because less magnetic free energy is available. Ni et al. present current sheets 
with lengths of order ~1 Mm which is a significant fraction of the size of the chromosphere, in 
contrast to lengths of ~10 km in our simulations. The thicknesses of the Ni et al. current sheets 
are generally much larger than A m which suggests that the neutrals and ions are strongly coupled 
in their regime. While recombination plays an important role in the simulations presented by L12, 
L13, and this work, the high temperatures found by Ni et al. make it less likely that net recombi- 
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Reconnection Rate 



Fig. 8.— The reconnection rate as a function of time as measured by the maximum change in 
magnetic flux among all X-points in the current sheet, &A z /&t. This rate is presented in dimen¬ 
sionless units according to Section 2.1; consequently, this rate has not been renormalized to values 
immediately upstream of the reconnection layer. 

nation will play an important role in the plasma continuity equation. Physical parameters such as 
magnetic field strength, temperature, density, and ionization fraction vary significantly within the 
chromosphere even at a single height, so it is likely that there are regions where each parameter 
regime is valid. These differences highlight the need for detailed parameter studies on how partial 
ionization impacts the reconnection process including the plasmoid instability at various locations 
in the chromosphere. 


3.4. Reconnection Rate 

The reconnection rate for all five simulations is shown in Fig. 8. This rate was measured as the 
maximum of the time derivatives of A z at all X-points within the reconnection region. A common 
alternative method for measuring the reconnection rate is to measure the inflow velocity divided by 
the outflow velocity; however, during asymmetric reconnection there is ambiguity in how to define 
the inflow velocity. The reconnection rate before t = 5 is not shown because an electric field was 
applied to initialize reconnection at early times. When interpreting this figure, it is important to 
recall that the initial total pressure and average upstream magnetic energy density are constant 
between runs. With this convention, magnetic and temperature asymmetry reduce the reconnection 
rate for all cases studied. 

From early times until around when plasmoids form, the reconnection rate is highest at the 
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X-point along x = 0. The reconnection rate increases as the current sheet thins. Around the time 
that plasmoids form, the reconnection rate at this X-point decreases. For some cases, this X-point 
bifurcates into an O-point. The peak reconnection rate then occurs at one of the nearby X-points 
that is not located along x = 0. Cases A, B, and D show a considerable increase in the reconnection 
rate after the formation of plasmoids. It is likely that Cases C and E would also show a similar 
increase in the reconnection rate; however, the simulations ended before this occurred. At even 
later times, structure on small enough scales could develop so that Hall reconnection could become 
important (see also Daughton et al. 2009; Shepherd & Cassak 2010). 


3.5. X-Point Motion and Flows Across the X-Point 

We consider three velocities that are relevant to the small-scale physics near the X-point located 
along the symmetry axis at x n = (0, y n ). The first velocity is the rate of motion of the X-point 
along the inflow direction: ^ L . This velocity is not a fluid velocity, but rather the velocity of a 
topologically stable magnetic feature (Murphy 2010). The second is the ion flow at the X-point 
along the inflow direction: Vi y (y n ). The third is the neutral flow at the X-point along the inflow 
direction: V ny (y n ). Differences between ^ and Vi y (y n ) must be due to non-ideal behavior such 
as resistivity or the Hall effect. Differences between Vi y (y n ) and V ny (y n ) correspond to momentum 
transfer between ions and neutrals (e.g., Eq. 9). 

Figure 9 shows Vi y {y n ), and V ny (y n ) as a function of time for all of the cases with asym¬ 
metric inflow. Case A is not shown because these velocities equal zero due to symmetry. The bump 
in velocity around t = 16 is a consequence of the initial conditions being slightly out of equilibrium; 
the resulting pulse is mostly but not entirely damped by a viscous layer at the boundary and the 
reflected pulse returns to the current sheet region at this time. The magnitude of this pulse is weak 
compared to the ion inflow speeds, but is the same order of magnitude as the velocities associated 
with each null point. For the rest of this section, we consider t > 20 for which reconnection is well 
established. 

In general, there is a small but nonzero difference between and Vi y (y n ). This similarity 
indicates that the magnetic field is predominantly carried by the ion flow. The residual difference 
in flow results from resistive diffusion of the magnetic field and the Hall effect. The difference 
between either of these velocities and the neutral flow at the X-point V ny (y n ) is greater, which 
is consistent with decoupling between the ion and neutral inflows. The V ny profile along the 
inflow direction remains roughly constant between Cases A and B which both have symmetric 
magnetic Helds but different temperature profiles, and also between Cases C-E which have the 
same asymmetric magnetic field configuration but differing initial temperature profiles. Case B has 
symmetric upstream magnetic fields but the y < 0 upstream region is warmer and less dense than 
the y > 0 upstream region. That « Vi y (y n ) < 0 indicates that the X-point is being carried 
primarily by the ions in the direction of the warmer upstream region. 
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Case B: (7.9, 7.9) G; (8750, 9250) K 



Case D: (10, 5) G; (8750, 9250) K 


Case C: (10, 5) G; (9000, 9000) K 



Case E: (10, 5) G; (9250, 8750) K 



Time 


Time 


Fig. 9.— A comparison between the X-point velocity < ^ i , the ion flow at the X-point along the 
inflow direction Vi y (y n ), and the neutral flow at the X-point V ny (y n ) for Cases B-E. We consider 
only the X-point along x = 0 and times after the early electric field application and before the 
onset of the plasmoid instability. Case A is not shown because all of these velocities equal zero due 
to symmetry. 

The simulations with magnetic asymmetry (Cases C, D, and E) show neutral flow from the 
weak field upstream region into the strong field upstream region. This flow results from a neutral 
pressure gradient that is pushing the neutrals from the weak field side into the strong field side 
where the neutral pressure is lower. This flow is able to occur because of imperfect coupling between 
species on scales below or comparable to the neutral-ion mean free path. Neutrals swept along with 
the outflow will have preferentially originated from the weak field upstream region. This flow may 
have important observational consequences because the abundances of neutrals swept along with 
the outflow are likely to better match the weak field upstream region. In contrast, there is net ion 











































flow inside the current sheet from the strong field upstream region into the weak field upstream 
region. While the neutrals are pushed by a neutral pressure gradient, the dominant force on the 
ions is a magnetic pressure gradient. 


4. DISCUSSION 

In this paper we perform simulations of asymmetric magnetic reconnection in partially ionized 
chromospheric plasmas using the plasma-neutral module of the HiFi framework. We include cases 
with symmetric or asymmetric upstream temperatures and magnetic field strengths. The plasma 
and neutrals are modeled as separate fluids with momentum and energy transfer between species 
due to collisions and charge exchange. These simulations self-consistently include ionization and 
recombination without assuming ionization equilibrium. 

Several of the properties of asymmetric partially ionized reconnection are qualitatively similar 
to the symmetric cases (Case A; L12; L13). The current sheet rapidly thins so that its thickness 
becomes comparable to the neutral-ion mean free path evaluated inside the current sheet. The 
ion and neutral outflows are strongly coupled, but the ion and neutral inflows are decoupled. 
The reconnecting magnetic field drags ions into the current sheet which leads to a significant 
enhancement of ion density so that the current sheet is out of ionization equilibrium. Recombination 
becomes of comparable importance to the outflow in the equation for ion conservation of mass. In 
this paper, we show that these properties are modified by the temperature and magnetic field 
asymmetries. 

During our simulations of asymmetric reconnection, the ion and neutral inflows remain decou¬ 
pled, but the decoupling is asymmetric. When there is magnetic asymmetry, there is net neutral 
flow through the current sheet from the weak magnetic field upstream region into the strong field 
upstream region. This neutral flow results from a large scale neutral pressure gradient and imper¬ 
fect coupling between ions and neutrals along the inflow direction. Similarly, the greater Lorentz 
force acting on the ions from the strong field upstream region leads to the ions pulling the X-point 
into the weak held upstream region. These effects are not present during symmetric simulations. 
An observational consequence of these neutral hows through the current sheet is that most of the 
neutrals swept along with the outflow are more likely to have originated from the weak magnetic 
held upstream region. This may be especially important if there are different elemental abundances 
in each upstream region (e.g., different FIP enhancements). 

Prior simulations of asymmetric reconnection in fully ionized plasmas have shown non-ideal 
hows through null points (e.g., Oka et al. 2008; Murphy 2010; Murphy et al. 2012). These hows 
result from non-ideal effects such as resistivity and the Hall effect. In resistive MHD, for example, 
the motion of null points results from a combination of resistive diffusion of the magnetic held and 
advection by the bulk plasma how (Murphy 2010). During asymmetric reconnection in partially 
ionized plasmas, an additional contributor to non-ideal plasma hows across null points is the electric 
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field contribution from electron-neutral drag (see the fourth term on the right hand side of Eq. 13). 
Neutral flow through the current sheet is a new effect that is not present in fully ionized situations. 

We include the Hall effect in our simulations to test whether or not a transition to Hall reconnec¬ 
tion can develop in this regime. We find that the out-of-plane quadrupole magnetic field associated 
with antiparallel Hall reconnection does develop. The out-of-plane magnetic field strength is a 
significant fraction of the in-plane field strengths just upstream of the current sheet. However, the 
characteristic low aspect ratio geometry and high reconnection rate expected during Hall recon¬ 
nection do not develop over the course of these simulations, including during the early nonlinear 
evolution of the plasmoid instability. The transition to Hall reconnection in a partially ionized 
plasma has been investigated analytically by Malyshkin & Zweibel (2011), and will continue to be 
investigated by ongoing simulation efforts to determine the conditions under which such a transition 
is likely to occur. 

After an initial phase, each of our simulations show the development of the plasmoid instability. 
These simulations capture the early nonlinear evolution until structures develop on scales compara¬ 
ble to the resolution scale. In cases with magnetic asymmetry, the plasnroids develop preferentially 
into the weak field upstream region, which is similar to fully ionized simulations. Recombination 
due to the enhanced ion density remains an important loss term in the plasma continuity equation. 
We anticipate that secondary merging of plasmoids will be further modified by these asymmetries. 
Our simulations complement the work of Ni et al. (2015) who use a single fluid framework to 
investigate the long term nonlinear evolution of the plasmoid instability in a different parameter 
regime. 

We investigate the coupling between magnetic and temperature asymmetries by comparing 
cases with the same magnetic asymmetry but different temperature asymmetries. We find that 
reconnection develops more quickly, occurs at a faster rate, and becomes unstable to plasmoid 
formation as we increase the temperature in the strong field upstream region while decreasing the 
temperature in the weak field upstream region. This change corresponds to Aincreasing in both 
upstream regions. The neutral flow through the current sheet does not noticeably change and 
therefore depends more strongly on the magnetic asymmetry than the temperature asymmetry. 

There are many important remaining directions for modeling work to understand how recon¬ 
nection occurs in the solar chromosphere and other weakly ionized plasmas. The simulations of L12, 
L13, and this paper have all simulated two-dimensional, antiparallel reconnection. The simulations 
presented in these works have modeled the early nonlinear evolution of the plasmoid instability 
in a weakly ionized plasma, and to model the later evolution will likely require a combination 
of increased resolution (including along the outflow direction to capture plasmoid merging) and 
increased diffusion. An alternative strategy would be to evolve the logarithm of density instead 
of the density itself (e.g., Lee et al. 2014) which precludes the possibility of the number density 
becoming negative. Ni et al. (2015) adopt adaptive mesh refinement to ensure adequate resolu¬ 
tion. The presence of a guide field can suppress the thinning of current sheets due to ambipolar 
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diffusion (Brandenburg &; Zweibel 1994, 1995), and thus should be considered in future efforts (see 
also Ni et al. 2015). Physical conditions in the chromosphere vary significantly even at a single 
height, so a detailed examination of parameter space will be necessary to characterize the different 
regimes of partially ionized reconnection in the chromosphere. The parameter studies may also 
include regimes relevant to partially ionized plasmas in the laboratory and in astrophysics. Future 
studies should investigate the impact of a realistic geometry and the interplay between small-scale 
physics and global dynamics (e.g., how small and large scales feed back on each other). Finally, 
much remains to be understood about the role of three-dimensional effects during weakly ionized 
reconnection and the behavior of current sheet thinning in cases with and without null points. 

In addition to the modeling efforts, much work remains to compare the simulation results to 
observations of chromospheric reconnection and to validate the results against laboratory experi¬ 
ments on partially ionized reconnection. A challenge in comparing simulations against observations 
is the disparity in length scales. The current sheets in these simulations have characteristic thick¬ 
nesses of < 10 2 m and lengths of < 10 km, while the spatial resolution of IRIS observations is of 
order 200 km, so direct diagnostics of the reconnection layer itself remain extremely challenging. 
However, the simulations could be used to predict velocities, spectra (including the charge state 
distributions of minor ions), densities, and morphologies that could then be compared against ob¬ 
servation. The experiments that have recently been performed at MRX provide an opportunity to 
directly validate HiFi’s plasma-neutral module against laboratory measurements where key quanti¬ 
ties can be observed using in situ probes. Using a code that has been validated against experimental 
data will improve confidence that it is able to accurately capture the relevant physical processes in 
the lower solar atmosphere and in astrophysical plasmas. 
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